      subroutine pd2dLiaoMcMechan(omega)


      use globais
      implicit none

      real*8 :: omega

      complex :: c1, c2, c3, c4, c5 ! variaveis auxiliares


      complex, dimension(:,:), allocatable :: LwB ! matriz do sistema em banda
!       para o solver do lapack
      integer :: NRHS, Neq, info
      integer :: KU,KL,LDAB,L
      integer, dimension(:), allocatable :: IPIV !auxiliar para o solver


      complex, dimension(:), allocatable :: f ! vetor do sistema  



      integer :: linha, col1,col2,col3,col4,col5 !  variáveis auxiliares para a montagem da matriz do sistema
      integer :: ndonoh, col !  variáveis auxiliares para prescrever superficie livre

      real*8:: xmin,zmin
      integer :: i,j
      call modelo2dcamadasviscoacoustico()

      allocate(creal(nx,nz))
      allocate(cimag(nx,nz))
      creal = real(c)
      cimag = aimag(c)

      xmin = -(Npml-1)*delta
      zmin = 0.00
      call systempsimage(  creal,nx,nz,delta,delta,zmin,xmin,"veloR",'velocidade (m/s)    ','profundidade (m)    ',&
		&'distancia (m)       ',wbox,Lz0/4.,Lx0/4.)

      call systempsimage(  cimag,nx,nz,delta,delta,zmin,xmin,"veloI",'velocidade (m/s)    ','profundidade (m)    ',&
		&'distancia (m)       ',wbox,Lz0/4.,Lx0/4.)

      call systempsimage(rho,nx,nz,delta,delta,zmin,xmin,"densi",'densidade  (kg/m3)  ','profundidade (m)    ',&
		&'distancia (m)       ',wbox,Lz0/4.,Lx0/4.)

      deallocate(creal)
      deallocate(cimag)



      Neq = Nx*Nz
      KL = Nx
      KU = Nx
      LDAB = 2*KL+KU+1
      allocate(LwB(LDAB,Neq))

      allocate(f(Neq))

!       "staggered-grid five-point stencil" 
! 	pags 1273: mixed-grid and staggered-grid finite-difference methods for 
!       frequency-domain acoustic wave modelling (hustedt, operto, virieux, 2004)




! 	write(*,*)'#Montando sistema...',kw,'/',Nomega
	do i=1, Nx
	    do j=1, Nz

		linha = (j-1)*Nx + i   !numero do noh ij

		col1 = (j-1)*Nx + i	   ! corresponde ao noh i,j  !numero do noh ij
		col2 = (j-1)*Nx + (i-1)  ! corresponde ao noh i-1,j !numero do noh da esquerda
		col3 = (j-1)*Nx + (i+1)  ! corresponde ao noh i+1,j	!numero do noh da direita
		col4 = ((j-1)-1)*Nx + i ! corresponde ao noh i,j+1  !numero do noh de cima
		col5 = ((j+1)-1)*Nx + i  ! corresponde ao noh i,j-1  !numero do noh de baixo
! 		write(*,*)'#   ',i,j,col1,col2,col3,col4,col5


		


		if(i.ne.1 .and. i.ne.Nx .and.j.ne.1 .and. j.ne.Nz )then

		  c2=1.
		  c3=1.
		  c4=1.
		  c5=1.
		  c1 = omega**2*delta**2/c(i,j)**2 - 2.

! 		  Lw(linha,col1) = c1
! 		  Lw(linha,col2) = c2
! 		  Lw(linha,col3) = c3
! 		  Lw(linha,col4) = c4
! 		  Lw(linha,col5) = c5
		  LwB(KL+KU+1+linha-col1,col1) = c1
		  LwB(KL+KU+1+linha-col2,col2) = c2
		  LwB(KL+KU+1+linha-col3,col3) = c3
		  LwB(KL+KU+1+linha-col4,col4) = c4
		  LwB(KL+KU+1+linha-col5,col5) = c5

		elseif(i.eq.1)then
		  c1 = -(1.0/delta + imag * omega/c(i,j))
		  c3 = +(1.0/delta)

		  c1 = 1.
		  c3 = 1.
! 		  Lw(linha,col1) = c1
! 		  Lw(linha,col3) = c3
		  LwB(KL+KU+1+linha-col1,col1) = c1
		  LwB(KL+KU+1+linha-col3,col3) = c3


		elseif(i.eq.Nx)then
		  c1 = +(1.0/delta + imag * omega/c(i,j))
		  c2 = -(1.0/delta)

		  c1 = 1.
		  c2 = 1.
! 		  Lw(linha,col1) = c1
! 		  Lw(linha,col2) = c2
		  LwB(KL+KU+1+linha-col1,col1) = c1
		  LwB(KL+KU+1+linha-col2,col2) = c2
		elseif(j.eq.1)then
		  c1 = -(1.0/delta + imag * omega/c(i,j))
		  c5 = +(1.0/delta)
		  c1 = 1.
		  c5 = 1.
! 		  Lw(linha,col1) = c1
! 		  Lw(linha,col5) = c5
		  LwB(KL+KU+1+linha-col1,col1) = c1
		  LwB(KL+KU+1+linha-col5,col5) = c5

		elseif(j.eq.Nz)then
		  c1 = +(1.0/delta + imag * omega/c(i,j))
		  c4 = -(1.0/delta)
		  c1 = 1.
		  c4 = 1.
! 		  Lw(linha,col1) = c1
! 		  Lw(linha,col4) = c4
		  LwB(KL+KU+1+linha-col1,col1) = c1
		  LwB(KL+KU+1+linha-col4,col4) = c4
		endif
      
		f(linha)=0.

	    enddo
	enddo



	! para prescrever u=0 na superficie livre (j=1), substituir a linha e a coluna referente ao noh ij inteira por zero, exceto o elemento da diagonal que é feito = a 1, zerar no vetor independente a linha de cada no da superficie livre  
	j=1 !nohs da superficie livre
	do i=1,Nx
	    ndonoh  = (j-1)*Nx + i    ! numero do noh i,j

	    linha = ndonoh
	    do col=1,Neq  ! zerando todos os elementos da linha "ndonoh"
		if((KL+KU+1+linha-col).ge.1 .and. (KL+KU+1+linha-col).le.LDAB .and. col .ge. 1 .and. col .le. Neq)then
		    LwB(KL+KU+1+linha-col,col)=0.
		endif
	    enddo

	    col=ndonoh
	    do linha=1,Neq  ! zerando todos os elementos da coluna "ndonoh"
		if((KL+KU+1+linha-col).ge.1 .and. (KL+KU+1+linha-col).le.LDAB .and. col .ge. 1 .and. col .le. Neq)then
		    LwB(KL+KU+1+linha-col,col)=0.
		endif	
	    enddo

	    linha = ndonoh	
	    col=ndonoh
	    LwB(KL+KU+1+linha-col,col)=1.    !coloca 1 no elemento da diagonal	    
	    f(linha)=0.                      !coloca 0 no vetor independente
	enddo
	!fim da condicao de contorno superficie livre





! 	open(1,file='matriz.txt')
! 	do i=1,Neq
! 	    write(1,2)(real(Lw(i,j)), j=1,Neq)
! 	enddo
! 	close(1)
!  2      format(9f10.4)


	f((jfonte -1)*Nx + ifonte) = Sw(kw)
! 	write(*,*) 'fonte',f((jfonte -1)*Nx + ifonte)

	! Solver: eliminacao de  gauss do lapack
	NRHS = 1
	allocate(IPIV(Neq))



! 	write(*,*)'#Resolvendo sistema........'
	call CGBSV(Neq,KL,KU,NRHS,LwB,LDAB,IPIV,f,Neq,INFO)! matriz em banda
! 	call CGESV(Neq,NRHS,Lw,Neq,IPIV,f,Neq,INFO)
	deallocate(IPIV)



 1      format(i10,2E20.6)

	do i=1,Nx
	    do j=1,Nz
		linha = (j-1)*Nx + i ! numero do noh i,j
		P(i,j,kw) =  f(linha)
	    enddo
	enddo


	do i=1,Nrec
! 		linhaReceiver = (jreceiver-1)*Nx + ireceiver
		linha = (ijReceiv(i,2)-1)*Nx + ijReceiv(i,1)  ! numero do noh do receptor i
		Preceiv(i,kw) =  f(linha)
	enddo		





      deallocate(LwB)
      deallocate(f)

      end subroutine pd2dLiaoMcMechan